
clear;
%insert "Data" file location in [...]
folder = '[...]Data';
cd(folder)
addpath(folder)

para0found = []
para1found = []

hrl0all = []
trl0all = []
hrl1all = []
trl1all = []

tch0 = []
tch1 = []
tnch0 = []
tnch1 = []

tchall0 = []
tchall1 = []

rhealthier0all = []
rhealthier1all = []
rtastier0all = []
rtastier1all = []

right0all=[]
right1all=[]


varRall=[]
Nall=[]

for i =  [13:25 27:29 31:32 34:37 39:55 57:60 62:65 66:68]
    X=readtable (strcat('Efood',  string(i), '_answers_new', '.xlsx'));
health=X.('Health')
taste=X.('Taste') 
nH=nansum(sign(abs(health)+1))
nT=nansum(sign(abs(taste)+1))
varR=((nH-1)/(nH+nT-2))*(nanstd(abs(health)))^2 + ((nT-1)/(nH+nT-2))*(nanstd(abs(taste)))^2 + 2*corr(abs(health),abs(taste),'rows','complete')*nanstd(abs(health))*nanstd(abs(taste))*((nT-1)/(nH+nT-2))*((nH-1)/(nH+nT-2))
varRall=[varRall;varR*(nH+nT)]
Nall=[Nall;nH+nT]
end

    avstdR = sqrt(sum(varRall)/sum(Nall))
   
probitf0ch = []
probitf1ch = []

tprobitf0ch = []
tprobitf1ch = []

taprobitf0ch = []
taprobitf1ch = []

hprobitf0ch = []
hprobitf1ch = []


for i = [13:25 27:29 31:32 34:37 39:55 57:60 62:65 66:68]
 
    
T0=readtable (strcat('behavior',  string(i), '_output_results0', '.xlsx'));
T1=readtable (strcat('behavior',  string(i), '_output_results1', '.xlsx'));

type0 = T0.('Type')
type1 = T1.('Type')

for r=1:90
   if type0(r,1)==1|type0(r,1)==3 
    chall0(r,1)=1
   else
    chall0(r,1)=0
   end 
end

for r=1:90
   if type1(r,1)==1|type1(r,1)==3 
    chall1(r,1)=1
   else
    chall1(r,1)=0
   end    
end

ch0 = chall0
nch0 = 1-chall0
ch0(chall0==0)=NaN
nch0(chall0==1)=NaN

ch1 = chall1
nch1 = 1-chall1
ch1(chall1==0)=NaN
nch1(chall1==1)=NaN

tch0 = [tch0, ch0]
tch1 = [tch1,ch1]
tnch0 = [tnch0, nch0]
tnch1 = [tnch1, nch1]

tchall0 = [tchall0, chall0]
tchall1 = [tchall1, chall1]

sel0=T0.('Selection') 
hrl0=(T0.('Stimuli_right_Health') - T0.('Stimuli_left_Health')) /avstdR
trl0=(T0.('Stimuli_right_Taste') - T0.('Stimuli_left_Taste')) /avstdR
sel1=T1.('Selection')
hrl1=(T1.('Stimuli_right_Health') - T1.('Stimuli_left_Health')) /avstdR
trl1=(T1.('Stimuli_right_Taste') - T1.('Stimuli_left_Taste')) /avstdR

for r=1:90
if sel0(r,1)=="right"
    right0(r,1)=1
elseif sel0(r,1)=="left"
    right0(r,1)=0
else
    right0(r,1)=NaN
end
end

for r=1:90
if sel1(r,1)=="right"
    right1(r,1)=1
elseif sel1(r,1)=="left"
    right1(r,1)=0
else
    right1(r,1)=NaN
end
end

rhealthier0 = sign(hrl0)
rhealthier1 = sign(hrl1)
rtastier0 = sign(trl0)
rtastier1 = sign(trl1)

hrl0all = [hrl0all,hrl0]
trl0all = [trl0all,trl0]

hrl1all = [hrl1all,hrl1]
trl1all = [trl1all,trl1]

rhealthier0all = [rhealthier0all,rhealthier0]
rhealthier1all = [rhealthier1all,rhealthier1]
rtastier0all = [rtastier0all, rtastier0]
rtastier1all = [rtastier1all, rtastier1]

right0all=[right0all,right0]
right1all=[right1all,right1]


mch0 = glmfit([hrl0 trl0 rhealthier0 rtastier0], right0, 'binomial','Link','probit')';
mch1 = glmfit([hrl1 trl1 rhealthier1 rtastier1], right1, 'binomial','Link','probit')';

probitf0ch = [probitf0ch;mch0]
probitf1ch = [probitf1ch;mch1]

if mch0(1,2)>0&mch0(1,3)>0
    mch0=mch0
elseif mch0(1,2)>0&mch0(1,3)<0
mch0 = glmfit([hrl0 rhealthier0 rtastier0], right0, 'binomial','Link','probit')'
mch0 = [mch0(1,1), mch0(1,2), 0, mch0(1,3), mch0(1,4)] 
elseif mch0(1,2)<0&mch0(1,3)>0
    mch0 = glmfit([trl0 rhealthier0 rtastier0], right0, 'binomial','Link','probit')'
    mch0 = [mch0(1,1), 0, mch0(1,2), mch0(1,3), mch0(1,4)] 
else
    mch0=[NaN,0,0,0,0]
end

if mch1(1,2)>0&mch1(1,3)>0
    mch1=mch1
elseif mch1(1,2)>0&mch1(1,3)<0
mch1 = glmfit([hrl1 rhealthier1 rtastier1], right1, 'binomial','Link','probit')'
mch1 = [mch1(1,1), mch1(1,2), 0, mch1(1,3), mch1(1,4)] 
elseif mch1(1,2)<0&mch1(1,3)>0
    mch1 = glmfit([trl1 rhealthier1 rtastier1], right1, 'binomial','Link','probit')'
    mch1 = [mch1(1,1), 0, mch1(1,2), mch1(1,3), mch1(1,4)] 
else
    mch1=[NaN,0,0,0,0]
end


%final probit coefficients (with positive Hrl and Trl coefficients)
tprobitf0 = [tprobitf0ch;mch0]
tprobitf1 = [tprobitf1ch;mch1]

alphaCh0 = mch0(1,2) / (mch0(1,2) + mch0(1,3))
sigmaCh0 = 1 / (mch0(1,2) + mch0(1,3))
biasCh0rl = mch0(1,1) / (mch0(1,2) + mch0(1,3))
biasCh0healthier = mch0(1,4) / (mch0(1,2) + mch0(1,3))
biasCh0tastier = mch0(1,5) / (mch0(1,2) + mch0(1,3))

alphaCh1 = mch1(1,2) / (mch1(1,2) + mch1(1,3))
sigmaCh1 = 1 / (mch1(1,2) + mch1(1,3))
biasCh1rl = mch1(1,1) / (mch1(1,2) + mch1(1,3))
biasCh1healthier = mch1(1,4) / (mch1(1,2) + mch1(1,3))
biasCh1tastier = mch1(1,5) / (mch1(1,2) + mch1(1,3))

para0found = [para0found;alphaCh0, biasCh0rl, sigmaCh0, biasCh0healthier, biasCh0tastier]
para1found = [para1found;alphaCh1, biasCh1rl, sigmaCh1, biasCh1healthier, biasCh1tastier]

hmch0 = glmfit([hrl0 rhealthier0 rtastier0], right0, 'binomial','Link','probit')';
hmch1 = glmfit([hrl1 rhealthier1 rtastier1], right1, 'binomial','Link','probit')';

hprobitf0ch = [hprobitf0ch;hmch0]
hprobitf1ch = [hprobitf1ch;hmch1]

tamch0 = glmfit([trl0 rhealthier0 rtastier0], right0, 'binomial','Link','probit')';
tamch1 = glmfit([trl1 rhealthier1 rtastier1], right1, 'binomial','Link','probit')';


taprobitf0ch = [taprobitf0ch;tamch0]
taprobitf1ch = [taprobitf1ch;tamch1]


end


%explore probit regressions

%first full model with biasRL and bias H/T

mch=[]

for o = 1:50
    
"-------SUBJECT---------two coeffs"  
o


"mch0"
mch0 =glmfit([hrl0all(:,o) trl0all(:,o) rhealthier0all(:,o) rtastier0all(:,o)], right0all(:,o), 'binomial','Link','probit')';
"mch1"
mch1 = glmfit([hrl1all(:,o) trl1all(:,o) rhealthier1all(:,o) rtastier1all(:,o)], right1all(:,o), 'binomial','Link','probit')';


"-------SUBJECT---------HD only"  
o


"mch0h"
mch0h =glmfit([hrl0all(:,o)  rhealthier0all(:,o) rtastier0all(:,o)], right0all(:,o), 'binomial','Link','probit')';
"mch1h"
mch1h = glmfit([hrl1all(:,o)  rhealthier1all(:,o) rtastier1all(:,o)], right1all(:,o), 'binomial','Link','probit')';


"-------SUBJECT---------TD only"  
o


"mch0t"
mch0t =glmfit([ trl0all(:,o) rhealthier0all(:,o) rtastier0all(:,o)], right0all(:,o), 'binomial','Link','probit')';
"mch1t"
mch1t = glmfit([ trl1all(:,o) rhealthier1all(:,o) rtastier1all(:,o)], right1all(:,o), 'binomial','Link','probit')';


%commands for checking the coeffs signs
%mch=[mch;mch0(1,2:3),mch1(1,2:3),mch0h(1,2),mch1h(1,2), mch0t(1,2),mch1t(1,2)]

end

%unidentifiable probit coeffs:
% 26 - cond 1; 21 - cond1

%coeffs positivity check (tables mch and mnch): all OK

%Use long dataset to run mixed probit (in Stata) to get group-level
%coefficients for the above subjects with unidentifiable probit coeffs for
%non-chall trials


%condition 0
 hrl_0 = .2447751  
 trl_0 = .224017    
 rhealthier_0 = .2350745 
 rtastier_0 = -.0262984
 const_0 = -.0411088    

alphanCh0_gr = (hrl_0/trl_0) / (1 + hrl_0/trl_0)
sigmanCh0_gr = alphanCh0_gr / hrl_0
biasnCh0rl_gr = const_0 * (alphanCh0_gr / hrl_0)
biasnCh0healthier_gr = rhealthier_0 * (alphanCh0_gr / hrl_0)
biasnCh0tastier_gr = rtastier_0 * (alphanCh0_gr / hrl_0)

%condition 1
 hrl_1 = .2160893
 trl_1 = .2037783  
 rhealthier_1 = .2683978 
 rtastier_1 = -.0310168
 const_1 = .0105055

 alphanCh1_gr = (hrl_1/trl_1) / (1 + hrl_1/trl_1)
sigmanCh1_gr = alphanCh1_gr / hrl_1
biasnCh1rl_gr = const_1 * (alphanCh1_gr / hrl_1)
biasnCh1healthier_gr = rhealthier_1 * (alphanCh1_gr / hrl_1)
biasnCh1tastier_gr = rtastier_1 * (alphanCh1_gr / hrl_1)

for o = [21 26]
   
    para0found(o,1) = alphanCh0_gr
    para0found(o,3) = sigmanCh0_gr
    para0found(o,2) = biasnCh0rl_gr
    para0found(o,4) = biasnCh0healthier_gr
    para0found(o,5) = biasnCh0tastier_gr
    
    para1found(o,1) = alphanCh1_gr
    para1found(o,3) = sigmanCh1_gr
    para1found(o,2) = biasnCh1rl_gr
    para1found(o,4) = biasnCh1healthier_gr 
    para1found(o,5) = biasnCh1tastier_gr
    
end


%macc is now for cut-off for 0
macc = 0.0000000000001


id=[]

prob_oo = []
probk=[]

    allsigmaH0=[]
    allsigmaT0=[]
    allsigmaH1=[]
    allsigmaT1=[]
    allro0=[ ]
    allro1=[]
      
    
for o = [1:50]

    id=[id;o]
    
Phj = []
Ptj=[]
scprobj=[] 
spprobj=[]

    jsigmaH0=[]
    jsigmaT0=[]
    jsigmaH1=[]
    jsigmaT1=[]
    jro0=[]
    jro1=[]
    
       
    for j = 0:0.001:1
           
     
     
    para0found_n(o,1) = j
    para0found_n(o,3) = macc
    
    para1found_n(o,1) = j
    para1found_n(o,3) = macc
    
       
%condition 0
if para0found(o,1) < macc & para0found_n(o,1) < macc
sigmaT0(o,1)=para0found(o,3)
sigmaT0_n(o,1)=para0found_n(o,3)
sigmaH0(o,1)= 10000
sigmaH0_n(o,1)= 10000
ro0(o,1)=0

elseif para0found(o,1) > 1 - macc & para0found_n(o,1) > 1 - macc
sigmaH0(o,1)=para0found(o,3)
sigmaH0_n(o,1)=para0found_n(o,3)
sigmaT0(o,1)= 10000
sigmaT0_n(o,1)= 10000
ro0(o,1)=0

elseif para0found(o,1) > 1 - macc & para0found_n(o,1) < macc
sigmaH0(o,1)=para0found(o,3)
sigmaH0_n(o,1)=10000
sigmaT0(o,1)= 10000
sigmaT0_n(o,1)= para0found_n(o,3)
ro0(o,1)=0

elseif para0found(o,1) < macc & para0found_n(o,1) > 1 - macc
sigmaH0(o,1)=10000
sigmaH0_n(o,1)=para0found_n(o,3)
sigmaT0(o,1)= para0found(o,3)
sigmaT0_n(o,1)= 10000
ro0(o,1)=0

else

%A means as in solution a) as in the text (Model section)

ro0A(o,1) = -sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1)))   
if ro0A(o,1)<0
ro0(o,1) = max(ro0A(o,1), 1/ro0A(o,1))
else
    ro0(o,1)=0
end
indA0(o,1) = 1 - abs( sign (ro0(o,1) - ro0A(o,1)) )  

sigmaH0f(o,1)=(para0found(o,3)/para0found(o,1))/sqrt(1-ro0(o,1)^2)
sigmaT0f(o,1)=(para0found(o,3)/(1-para0found(o,1)))/sqrt(1-ro0(o,1)^2)
sigmaH0f_n(o,1)=(para0found_n(o,3)/para0found_n(o,1))/sqrt(1-ro0(o,1)^2)
sigmaT0f_n(o,1)=(para0found_n(o,3)/(1-para0found_n(o,1)))/sqrt(1-ro0(o,1)^2)

if sigmaH0f(o,1)==inf
    sigmaH0f(o,1)=10000
else
sigmaH0f(o,1) = sigmaH0f(o,1)
end

if sigmaT0f(o,1)==inf
    sigmaT0f(o,1)=10000
else
sigmaT0f(o,1) = sigmaT0f(o,1)
end

if sigmaH0f_n(o,1)==inf
    sigmaH0f_n(o,1)=10000
else
sigmaH0f_n(o,1) = sigmaH0f_n(o,1)
end

if sigmaT0f_n(o,1)==inf
    sigmaT0f_n(o,1)=10000
else
sigmaT0f_n(o,1) = sigmaT0f_n(o,1)
end

if para0found(o,1) > 1 - macc
   fsigmaH0f(o,1) = macc
   fsigmaT0f(o,1) = 10000
elseif para0found(o,1) < macc
     fsigmaH0f(o,1) = 10000
     fsigmaT0f(o,1) = macc
else
fsigmaT0f(o,1)=sigmaH0f(o,1) * (abs(ro0(o,1)) / ((1-para0found(o,1))/para0found(o,1)))
fsigmaH0f(o,1)=sigmaT0f(o,1) * (abs(ro0(o,1)) / (para0found(o,1)/(1-para0found(o,1))))
end

if para0found_n(o,1) > 1 - macc
   fsigmaH0f_n(o,1) = macc
   fsigmaT0f_n(o,1) = 10000
elseif para0found_n(o,1) < macc
     fsigmaH0f_n(o,1) = 10000
     fsigmaT0f_n(o,1) = macc
else
fsigmaT0f_n(o,1)=sigmaH0f_n(o,1) * (abs(ro0(o,1)) / ((1-para0found_n(o,1))/para0found_n(o,1)))
fsigmaH0f_n(o,1)=sigmaT0f_n(o,1) * (abs(ro0(o,1)) / (para0found_n(o,1)/(1-para0found_n(o,1))))
end

sigmaH0(o,1)   = indA0(o,1)*fsigmaH0f(o,1)  + (1-indA0(o,1))*sigmaH0f(o,1)         %(1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * fsigmaH0f(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * sigmaH0f(o,1)
sigmaH0_n(o,1) = indA0(o,1)*sigmaH0f_n(o,1) + (1-indA0(o,1))*fsigmaH0f_n(o,1)       %(1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * sigmaH0f_n(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * fsigmaH0f_n(o,1)
sigmaT0(o,1)   = indA0(o,1)*sigmaT0f(o,1)   + (1-indA0(o,1))*fsigmaT0f(o,1)         % (1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * sigmaT0f(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * fsigmaT0f(o,1)
sigmaT0_n(o,1) = indA0(o,1)*fsigmaT0f_n(o,1)  + (1-indA0(o,1))*sigmaT0f_n(o,1)          % (1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * fsigmaT0f_n(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * sigmaT0f_n(o,1)

if sigmaH0(o,1)==0
     sigmaH0(o,1)=macc
else
sigmaH0(o,1)=sigmaH0(o,1)
end

if sigmaH0_n(o,1)==0
     sigmaH0_n(o,1)=macc
else
sigmaH0_n(o,1)=sigmaH0_n(o,1)
end

if sigmaT0(o,1)==0
     sigmaT0(o,1)=macc
else
sigmaT0(o,1)=sigmaT0(o,1)
end

if sigmaT0_n(o,1)==0
     sigmaT0_n(o,1)=macc
else
sigmaT0_n(o,1)=sigmaT0_n(o,1)
end



end

%condition 1
if para1found(o,1) < macc & para1found_n(o,1) < macc
sigmaT1(o,1)=para1found(o,3)
sigmaT1_n(o,1)=para1found_n(o,3)
sigmaH1(o,1)= 10000
sigmaH1_n(o,1)= 10000
ro1(o,1)=0

elseif para1found(o,1) > 1 - macc & para1found_n(o,1) > 1 - macc
sigmaH1(o,1)=para1found(o,3)
sigmaH1_n(o,1)=para1found_n(o,3)
sigmaT1(o,1)= 10000
sigmaT1_n(o,1)= 10000
ro1(o,1)=0

elseif para1found(o,1) > 1 - macc & para1found_n(o,1) < macc
sigmaH1(o,1)=para1found(o,3)
sigmaH1_n(o,1)=10000
sigmaT1(o,1)= 10000
sigmaT1_n(o,1)= para1found_n(o,3)
ro1(o,1)=0

elseif para1found(o,1) < macc & para1found_n(o,1) > 1 - macc
sigmaH1(o,1)=10000
sigmaH1_n(o,1)=para1found_n(o,3)
sigmaT1(o,1)= para1found(o,3)
sigmaT1_n(o,1)= 10000
ro1(o,1)=0

else
    
%ro1A(o,1) = -sqrt( (para1found_n(o,1)/(1-para1found_n(o,1))) * ((1-para1found(o,1))/para1found(o,1)))   
ro1A(o,1) = -sqrt( (para1found(o,1)/(1-para1found(o,1))) * ((1-para1found_n(o,1))/para1found_n(o,1)))   
if ro1A(o,1)<0
ro1(o,1) = max(ro1A(o,1), 1/ro1A(o,1))
else
    ro1(o,1)=0
end
indA1(o,1) = 1 - abs( sign (ro1(o,1) - ro1A(o,1)) ) 

sigmaH1f(o,1)=(para1found(o,3)/para1found(o,1))/sqrt(1-ro1(o,1)^2)
sigmaT1f(o,1)=(para1found(o,3)/(1-para1found(o,1)))/sqrt(1-ro1(o,1)^2)
sigmaH1f_n(o,1)=(para1found_n(o,3)/para1found_n(o,1))/sqrt(1-ro1(o,1)^2)
sigmaT1f_n(o,1)=(para1found_n(o,3)/(1-para1found_n(o,1)))/sqrt(1-ro1(o,1)^2)

if sigmaH1f(o,1)==inf
    sigmaH1f(o,1)=10000
else
sigmaH1f(o,1) = sigmaH1f(o,1)
end

if sigmaT1f(o,1)==inf
    sigmaT1f(o,1)=10000
else
sigmaT1f(o,1) = sigmaT1f(o,1)
end

if sigmaH1f_n(o,1)==inf
    sigmaH1f_n(o,1)=10000
else
sigmaH1f_n(o,1) = sigmaH1f_n(o,1)
end

if sigmaT1f_n(o,1)==inf
    sigmaT1f_n(o,1)=10000
else
sigmaT1f_n(o,1) = sigmaT1f_n(o,1)
end

if para1found(o,1) > 1 - macc
   fsigmaH1f(o,1) = macc
   fsigmaT1f(o,1) = 10000
elseif para1found(o,1) < macc
     fsigmaH1f(o,1) = 10000
     fsigmaT1f(o,1) = macc
else
fsigmaT1f(o,1)=sigmaH1f(o,1) * (abs(ro1(o,1)) / ((1-para1found(o,1))/para1found(o,1)))
fsigmaH1f(o,1)=sigmaT1f(o,1) * (abs(ro1(o,1)) / (para1found(o,1)/(1-para1found(o,1))))
end

if para1found_n(o,1) > 1 - macc
   fsigmaH1f_n(o,1) = macc
   fsigmaT1f_n(o,1) = 10000
elseif para1found_n(o,1) < macc
     fsigmaH1f_n(o,1) = 10000
     fsigmaT1f_n(o,1) = macc
else
fsigmaT1f_n(o,1)=sigmaH1f_n(o,1) * (abs(ro1(o,1))/ ((1-para1found_n(o,1))/para1found_n(o,1)))
fsigmaH1f_n(o,1)=sigmaT1f_n(o,1) * (abs(ro1(o,1)) / (para1found_n(o,1)/(1-para1found_n(o,1))))
end


sigmaH1(o,1)   = indA1(o,1)*fsigmaH1f(o,1)  + (1-indA1(o,1))*sigmaH1f(o,1)         %(1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * fsigmaH0f(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * sigmaH0f(o,1)
sigmaH1_n(o,1) = indA1(o,1)*sigmaH1f_n(o,1) + (1-indA1(o,1))*fsigmaH1f_n(o,1)       %(1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * sigmaH0f_n(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * fsigmaH0f_n(o,1)
sigmaT1(o,1)   = indA1(o,1)*sigmaT1f(o,1)   + (1-indA1(o,1))*fsigmaT1f(o,1)         % (1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * sigmaT0f(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * fsigmaT0f(o,1)
sigmaT1_n(o,1) = indA1(o,1)*fsigmaT1f_n(o,1)  + (1-indA1(o,1))*sigmaT1f_n(o,1)          % (1-abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) )) ) * fsigmaT0f_n(o,1)   + abs(sign(abs(ro0(o,1))-sqrt( (para0found(o,1)/(1-para0found(o,1))) * ((1-para0found_n(o,1))/para0found_n(o,1))) ))  * sigmaT0f_n(o,1)

if sigmaH1(o,1)==0
     sigmaH1(o,1)=macc
else
sigmaH1(o,1)=sigmaH1(o,1)
end

if sigmaH1_n(o,1)==0
     sigmaH1_n(o,1)=macc
else
sigmaH1_n(o,1)=sigmaH1_n(o,1)
end

if sigmaT1(o,1)==0
     sigmaT1(o,1)=macc
else
sigmaT1(o,1)=sigmaT1(o,1)
end

if sigmaT1_n(o,1)==0
     sigmaT1_n(o,1)=macc
else
sigmaT1_n(o,1)=sigmaT1_n(o,1)
end



end



       
    Ph0_o = []  
    Pt0_o = [] 
    Ph1_o = []  
    Pt1_o = []
    spprob0_o = []
    scprob0_o = []
    spprob1_o = []
    scprob1_o = []
         
        
     Ph0_ou = []  
    Pt0_ou = [] 
    Ph1_ou = []  
    Pt1_ou = []
    spprob0_ou = []
    scprob0_ou = []
    spprob1_ou = []
    scprob1_ou = []
    
     for i = [1:90]
       
         % differences for chosen minus unchosen item
        mu0 = [right0all(i,o)*hrl0all(i,o) + (1-right0all(i,o))*(-hrl0all(i,o)),right0all(i,o)*trl0all(i,o) + (1-right0all(i,o))*(-trl0all(i,o))]
        mu1 = [right1all(i,o)*hrl1all(i,o) + (1-right1all(i,o))*(-hrl1all(i,o)),right1all(i,o)*trl1all(i,o) + (1-right1all(i,o))*(-trl1all(i,o))]
                       
        
        sigma0 =  [sigmaH0(o,1), ro0(o,1); ro0(o,1),sigmaT0(o,1)]
        sigma1 =  [sigmaH1(o,1), ro1(o,1); ro1(o,1),sigmaT1(o,1)]
        
              
        if sigma0(1,1)>0&sigma0(1,2)>-1.01&sigma0(2,1)>-1.01&sigma0(2,2)>0&sigma0(1,1)<inf&sigma0(2,2)<inf
        Ph0_o = [Ph0_o; mvncdf([0 -inf],[inf inf],mu0,sigma0)]
        Pt0_o = [Pt0_o; mvncdf([-inf 0],[inf inf],mu0,sigma0)]
        spprob0_o = [spprob0_o; mvncdf([0 0],[inf inf],mu0,sigma0) / (mvncdf([0 0],[inf inf],mu0,sigma0) + mvncdf([-inf -inf],[0 0],mu0,sigma0))]
        scprob0_o = [scprob0_o; mvncdf([0 -inf],[inf 0],mu0,sigma0) / (mvncdf([0 -inf],[inf 0],mu0,sigma0) + mvncdf([-inf 0],[0 inf],mu0,sigma0))]
                  
                  
                                 
        else
            Ph0_o = [Ph0_o; NaN]
            Pt0_o = [Pt0_o; NaN]
            spprob0_o = [spprob0_o; NaN]
            scprob0_o = [scprob0_o; NaN]
          
        end
        
        
        if sigma1(1,1)>0&sigma1(1,2)>-1.01&sigma1(2,1)>-1.01&sigma1(2,2)>0&sigma1(1,1)<inf&sigma1(2,2)<inf
         Ph1_o = [Ph1_o; mvncdf([0 -inf],[inf inf],mu1,sigma1)]
        Pt1_o =  [Pt1_o; mvncdf([-inf 0],[inf inf],mu1,sigma1)]
        spprob1_o = [spprob1_o; mvncdf([0 0],[inf inf],mu1,sigma1) / (mvncdf([0 0],[inf inf],mu1,sigma1) + mvncdf([-inf -inf],[0 0],mu1,sigma1))]
        scprob1_o = [scprob1_o; mvncdf([0 -inf],[inf 0],mu1,sigma1) / (mvncdf([0 -inf],[inf 0],mu1,sigma1) + mvncdf([-inf 0],[0 inf],mu1,sigma1))]
               
       
                               
            
        else
            Ph1_o = [Ph1_o; NaN]
            Pt1_o = [Pt1_o; NaN]
            spprob1_o = [spprob1_o; NaN]
            scprob1_o = [scprob1_o; NaN]
            
        end
        
       
   
        
        
     end
     
     
   
     
          Ph = [Ph0_o; Ph1_o]
          Pt = [Pt0_o; Pt1_o]
          spprob = [spprob0_o; spprob1_o]
          scprob = [scprob0_o; scprob1_o]
          
          
   chall = [tchall0(:,o);tchall1(:,o)]
   
            
   
   
   
   cond = [zeros(90,1); ones(90,1)]
   prob_o = [o.*ones(180,1), cond, chall] 
   
Phj = [Phj, Ph]
Ptj = [Ptj, Pt]
scprobj = [scprobj, scprob]
spprobj = [spprobj, spprob]


    jsigmaH0=[jsigmaH0, sigmaH0(o,1)]
    jsigmaT0=[jsigmaT0, sigmaT0(o,1)]
    jsigmaH1=[jsigmaH1, sigmaH1(o,1)]
    jsigmaT1=[jsigmaT1, sigmaT1(o,1)]
    jro0=[jro0, ro0(o,1)]
    jro1=[jro1, ro1(o,1)]

  
    end  
    probk=[probk;nanmean(Phj,2), nanmean(Ptj,2),nanmean(scprobj,2),nanmean(spprobj,2)]
    prob_oo = [prob_oo; prob_o]
    
    
    allsigmaH0=[allsigmaH0; jsigmaH0]
    allsigmaT0=[allsigmaT0; jsigmaT0]
    allsigmaH1=[allsigmaH1; jsigmaH1]
    allsigmaT1=[allsigmaT1; jsigmaT1]
    allro0=[allro0;jro0 ]
    allro1=[allro1; jro1]
    
         
end

tprob = ["subject", "condition", "chall", "Ph","Pt", "scprob", "spprob"; prob_oo, probk]  
 
 
   msigmaH0=mean(allsigmaH0,2)
    msigmaT0=mean(allsigmaT0,2)
    msigmaH1=mean(allsigmaH1,2)
    msigmaT1=mean(allsigmaT1,2)
    mro0=mean(allro0,2)
    mro1=mean(allro1,2)

  
   
Parameters = ["alpha0", "ro0", "sigmaH0", "sigmaT0", "sigma0","bias0rl", "bias0healthier","bias0tastier", "alpha1","ro1", "sigmaH1", "sigmaT1", "sigma1", "bias1rl", "bias1healthier", "bias1tastier",; para0found(:,1), mro0, msigmaH0, msigmaT0, para0found(:,3),para0found(:,2), para0found(:,4),para0found(:,5),para1found(:,1), mro1, msigmaH1, msigmaT1, para1found(:,3),para1found(:,2), para1found(:,4),para1found(:,5) ]

save('Approach3_1000')